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1.0 INTRODUCTION 


Numerical stability is an important factor in evaluating the merit of a particu- 
lar numerical solution. It appears that users of numerical algorithms often as- 
sume that if their system of linear, constant coefficient differential equations 
has eigenvalues with real nonpositive parts, and since the (exact) solutions of 
initial value problems are stable, it is safe to assume that the numerical solu- 
tions are also stable. However, an increasing number of dynamical problems with 
eigenvalues (multiplicity of two or higher) situated close to the imaginary axis 
seem to present serious difficulties in regards to their numerical integration. 
For instance, in the study of orbital mechanics (e.g., the linearized 
Kustaanheimo-Stiefel elements equations for the unperturbed two-body problem) , 
the differential systems often have imaginary or zero eigenvalues of high 
multiplicity. Certainly a thorough understanding of the numerical solutions of 
the linear systems would be required prior to undertaking the difficult study of 
the stability of numerical solutions of nonlinear systems. 

While various texts and papers have addressed the stability regions for Runge- 
Kutta codes and stable solutions, it seems that the analysis is inadequate for 
many of the current problems. The heavy dependence on test equations that are 
scalar or (at best) diagonal systems is a part of the problem. In this paper 
the emphasis is on systems (excluding two preliminary sections) and 
illustrations. Theorems and definitions are directed towards systems with 
eigenvalues 


X = p (cos 6 + i sin 0) 


where P > o and 90° < 0 < 270° 
and i = /-T 


Concepts such as amplification matrices, numerical kernels, stable, and exponen- 
tially stable numerical solutions are examined. In section 5.0 the techniques 
of previous sections are applied to certain systems that have Jordan forms, 
which are nondiagonal with particular interest in the case of imaginary or zero 
eigenvalues. 


2.0 AMPLIFICATION FACTORS 


The general four-stage explicit Runge-Kutta method advances one step according 
to 


4 

y n+ 1 : Vn + h n ^ “i fi (1) 

i = 1 
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where 


i-1 

fi = f(t n + aih n , v n ♦ h n E b iJ*V ^ 

j=1 


for 1 = 1,2, 3, 4 and h n s t n+ i - t n . The usual approach of Batching Taylor se- 
ries expansions, Lambert (ref. 1) implies that 


i-1 

ai = o, ai = £ bij, i = 2,3,** (3) 

j=1 


and consistency requires that 


4 

£ = 1 ( ^ ) 
i= 1 


Numerical stability for Runge-Kutta schemes considers error propagation for the 
linear, scalar trial equation 

£ = \y (5) 

where X can be complex valued. Since f(t,y) = Xy, substitution in equation 

(2) gives 

fi = ^Yn 

f 2 s Xy n (l ♦ Xh n a2) 

f 3 = ^yn ( 1 + (Xh n ) + a2b32 (Xh n ) 2 ) 

f4 = Xy n [ 1 ♦ a4 (Xh n ) + (a2b42 ♦ a3b43)(Xh n )2 + a 2 b 32 b ^3 (Xh n )3) (6) 

where the equation (3) relationship is used. Consequently, evaluation of 
equation (1) using equation (6) yields 


y n +1 = P(Xh n ) y n 


(7) 
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where 


p(Xhn) * 1 


+ Ah n ♦ 




(Ah n ) 2 ♦ (a2b32d)3 + a2b42u>H ♦ a 3bj43u>i4 ) ( Xb n ) 3 


+ (a2b32b43a)4)(Ah n ) i, 


( 8 ) 


Equation (4) was used to determine the coefficient of xh n in equation (8). 

The scalar polynomial (eq. (8)) is called the amplification factor for the 
four-stage Runge-Kutta method. If the method is to have fourth order, then 

y(t n +i) - y n +i = 0(h5) (9) 

where y(t n+i ) is the exact solution at t n+i and y n+1 is the (n+i) th iter- 
ate of the algorithm. Solving equation (5) implies the exact relationship: 


Ah n 

y(t n+ i) * e y(t n ) 


( 10 ) 


where 


= -n+1 * ^n 

Subtracting equation (7) from equation (10) implies that 


y(tn+i) “ yn+i 


Ah n 

e 


y(t n ) - p(\h n ) y n 


( 11 ) 


But 


y(t n ) s y n ♦ 0(h5) 


and by substitution to 0(h5), 


y (^n+l ) ~ yn+1 



- P(Xh n )) y n 


( 12 ) 


3 
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Since to 0(h5) 


e ^ h n _ j 

i=o ** 


it can be deducted from equation (12) that the error amplification factor is 


p(Xh n ) = Z 
i=o 


(Xhn ) 1 

i! 


(13) 


provided the Runge-Kutta scheme is to be fourth order. The region of absolute 
stability in this particular method is the set of all Xh for which 


i p(Xh) ( < 1 


( 1*0 


where h is a positive number. That is, the graph of the region (in the com- 
plex plane) 


S 


[u : u complex and 



< 


(15) 


is the stability region of all fourth-order explicit Runge-Kutta methods (with 
four stages). Figure 1 shows the graph of S. 

For general Runge-Kutta methods, the amplification factor p(y) is a polynomial 
in y = Xh, which is a function of the order q and the number of stages s of 
the method 


q p i s 

p(p) = Z — ♦ Z Hi U 1 

i =0 11 i=q+1 


(16) 


where the are functions of the Runge-Kutta parameters. Additional informa- 

tion for s > q > *4 is given in references 2 (sec. 3.3) and 3. 


Suppose that 
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y( fc n) " - e n ^7) 

where e n denotes the global error at tine t n , then from equation (11) 

Xh n 

y(trn-i) - y n +i = • y(t n > ■ p (Xh n> yn 


X h 

£ n +l * « n y(t n > - P(^h n ) (y(t n ) - e n ) » 

e n+1 = P( x ^n^ e n ^®) 

to CKh^). Equation (18) Justifies the term amplification factor, and it is 
clear that if | p(Xh n ) j < 1 , then one-step error control has been achieved. 

Remarks. The ooncept of a numerical stability region is dependent on (a) the 
method (e.g., Runge-Kutta or linear multistep, implicit or explicit, order, 
number of stages, etc.), and (b) the nature of the trial equation. Regarding (b), 
it is important to realize that while the region is independent of eigenvalues 
per se, the analysis that leads to the amplification factor was very dependent 
on the type of equation (i.e., eq. (5)). 

In the "positive lobes" of the stability region S, the set 
LcS, L = (yeS:y complex arid Rey > o ) 


where Rey is the real part of the complex number y, is an anomaly associated 
with certain Runge-Kutta schemes. This phenomenon occurs here because in the 
partial series expansion of 


Xh ah+iBh 

e : e 


where 


i^ = -1 and X = a + iB 


The condition 


5 
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4 

Z 

J=o 


(a ♦ i0)J 

Jj 


I * 1 


can occur (because of crossover terms) with a > o where 


cxh+i8h 


1 


if and only if a=o (assuming h>o) . The L is usually deleted from the graph 
of the stability region S for the following reason: pick any y with 
Rey < o and consider the one-half line or ray starting at the origin and 
passing through y. This ray is divided into two sets, each of which is connected 
namely, points of the ray inside S and those exterior to S. How consider a 
y (near the imaginary axis with Rey > o) such that its ray passes through L. 

The set L partitions this ray into two sets but the points of the ray in the 

unstable regions are disconnected, and in fact, separated by the set L. This 
unstable region near the origin is what detracts from the applicability of 
L as part of the stability region. See figure 7 for the plot of | p(Xh) j 

with Xi = Pi (cos 85° + i sin 85°) and p-| = 2, P 2 = 5, and P 3 = 10. Figure 

7 illustrates the nature of L. 


3-0 SYSTEMS OF EQUATIONS 

Equations (1) and (2) of the Runge-Kutta method remain unchanged for systems ex- 
cept that y and f are to be considered as n dimensional vectors. The 
Runge-Kutta parameters and the step sizes are scalars. The linear trial system 
is 

l = i|> l (19) 


where is an n by n constant, diagonal matrix 


= diag { X 1 ,X 2 , . • • ,X n } ( 20 ) 

where ^ = (yi,.*.,y n ) T and Xi are possibly multiple. 

Proceeding as before, the relationship 


Zrn+1 r 1m 


(21) 


6 
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is obtained where P is an n by n diagonal matrix, 


P(Xh B ) * diag { p(X 1 h B ), p(X 2 h m ) , . . . ,p(X n h B )> (22) 


The scalar p(Xih B ) is defined by equation (8) for each i. The matrix 
P(Xh m ) = PChjj^) will be referred to as the amplification matrix. The global 
error (eq. (17)) and error propagation formula (eq. 08)) are as before, 
except gjn is now an n component vector. The region of absolute stability 
is unchanged; however, to conclude that the numerical integration proceeds 
in a stable way, it will be necessary that h n be selected sufficiently 
small so that 


I p(Xjh m ) j 5 . ^ > X - 1 , . • * ,n 


4.0 NUMERICAL STABILITY 

The concept of distance will be significant in this section. The norm |j ||, 
defi ned 


li Ax II 

- sup ■ — ■ 

2L * 9 . || £ ! | 


for square matrices A, and either 


|| x |j = max Ix^ or ||x|| = ( Z Xi 2 ) 1/2 
i=1,...,n 1=1 


for vectors will be adopted. If A is any square matrix, then it follows from 
the definition of || A || that 

11 Ax 1| < | | A| J | |x| | , all x 


Numerical stability is a hybrid concept in that it is both method and equation 
dependent. In contrast to Lyapunov stability, it would seem unreasonable to ex- 
pect that a classification of eigenvalues of a linear system would suffice as 
characteristic (of numerical sta - .tty) although such a possibility is not 
ruled out. What is wanted is a definition of stability for a numerical solution 


II * || 


7 
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, (hi)”" 1 ) 
iso iso 


which has resulted from the implementation of a numerical algorithm for an 
initial value problem 


i = £(t,z) I 

£(t 0 ) = loi < 2 3) 


The definition should enable the user to classify the given solution on the 
basis of the behavior of certain nearby numerical solutions using the same 
steps and the same algorithm. The need to keep the steps invariant is required 
to avoid confusing stability and convergence. Tne concept of numerical stabil- 
ity introduced below is global (from a numerical viewpoint); i.e., based on 
behavior at discrete points in the interval (t 0 ,tf) = { t : t 0 < t < tf) where 


N-1 

tf = t 0 + I hi 
iso 


and makes use of a "discrete tube." Two basic assumptions are made; first, 
equation (23) has unique solutions in the region of interest and secondly, 
roundoff error does not occur (i.e., infinite machine precision). 

Definition 1. A numerical solution 



to an initial value problem (eq. (23)) that is generated by a specific nurierical 
algoriihm is exponentially stable if positive numbers p and 6, o < p < 1 
exist , such that 


1 1 Hi - 5.i 1 1 < P 1 li io - 5.0 H 


(25) 


for all numerical solutions 


8 
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, N N-1 

i 2 i) » {hj.} 

i=o i=o 

with 

II Zo ~ 2 o II £ 5 (26) 

The numerical solution, equation (24), is stable if p = 1; i.e., 

I I Zi ” 2i I I £ I i Zo ” £o 1 1 » » * * • (27) 

N 

The solution {zj_} is computed according to the same algorithm as equation 
i=o 

(24). Equation (24) will be exponentially unstable if in every <5 neighborhood 

N 

£ 0 (eq. (26)) there exists a numerical solution {z^} and p>1 such that 

i=o 


II Zi - z i II > P 1 II Zo - zo II » i=1 , . . . ,N 


( 28 ) 


Discussion . The above categories of numerical stability are not mutually exclu- 
sive, (i.e.. an exponentially stable solution is stable but not necessarily con- 
versely, and there exist numerical solutions that are neither stable nor expon- 
entially unstable) . It will often be beneficial to view an initial value £ 0 

N-1 

and step sequence {h^} as given; the implementation of the algorithm 

i=o 

N 

results in { y-; } 


Geometrically the concepts of exponentially stable and unstable are straight- 
forward. Figure 2 illustrates the geometry. A discrete tube is a discon- 
nected set of N+1 disks in tx R n that are norm dependent for their geometry 

i-1 

and that are centered on at time tj_ = t 0 + Z hj. The radii of the 

j=o 

disks are identified with the right sides of equations (25) and (27). 

If the numerical solution, equation (24), is exponentially stable, then applica- 
tion of the given algorithm to any initial value z 0 within 6 distance of 


9 
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n 

y 0 will produoe a aequenae of veotors J ^ ^ auoh that z± ia within 

P 1 II lo " £o II distanoe of y ± ; for eaoh i, i*1,... f N. Figure 2 illuatratea 
thla aituation if | 1 - Zq || = 6 z Q actually falla on the circumfer- 

ence of the initial diak) then z^ would be in or on the i th disk and p<1. 

For atability, all the diaka have radii H ^ - Zq ||; thua z A deviates from 
li (for all i) by no more than z D from ^ 0 . 

Theorem 1 . Consider the initial value problem 


i = ^ . y(t Q ) = io 


(29) 


where = diag {X 1 ,... f X n ) with ReXj^ < o , for all i. 


In addition, assume that 




(24) 


is a numerical solution of equation (29) generated by some explicit Runge-Kutta 
algorithm with stability region S, 

S = {p:u complex and |p(p)| < 1} (30) 

where p is the amplification factor. Then there exists a number h*>o such 
that the numerical solution (eq. (24)) is stable if hie(o,h*} for all i. 

Proof . From equation (21), 


Kra+1 = P(Xh m ) 


(3D 


where 


P(Xhnj) = diag (p(X 1 h m ),...,p(X n h m )} 


10 
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for all m, m=o, 1 , . . . ,N-1 . The symbol PUhg,) is written as P ( h m 4> ) where 
the p(X i h m ) are defined in equation (8) for the general, fourth-order 
Runge-Kutta method and by similar amplification factors in general. Repeated 
iterations of equation (3D with m=0,1,2,... implies that for any i, 
i=o, 1 , . . . ,N-1 : 

Hi = 40 P(hi_ 2 40 . . .P(h 0 40 ) n 0 (32) 


Equation (32), which shows how a solution propagates, can be written (using the it 
product) as 


i-1 


Hi = l it P(hj 4')/ Ho 

\j=0 


(33) 


and as a consequence of the diagonal nature of the amplification matrices, 


i-1 


|i-1 i-1 ) 

n P(h} 40 = diag \ it p( X i h j ) , . . . , tr p(X n h j ) ^ 

j=0 (j=o j=o I 


(3D 


Similarly, for a numerical solution with initial value z^, it is concluded 
that 


Si 


0-1 1 

I n P(h 1 40 

\J=o / 


So 


(35) 


The matrix in equation (35) (i.e., the tt product) will be called the numerical 
kernel. Subtracting equation (35) from equation (33) and taking norms implies 
that 


i-1 


I I Hi " Si li = II ( 11 ^Ho “ Sc 

J=0 


(36) 


The k^ component of the vector 


^ tt p(hj 'I')! - z 0 ) equals ^ tt p( h jX k [ 


yko “ z ko 


11 
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where yj< 0 and z ko are the k^ 1 components of jTq and z 0 . Thus, 
ll^ 1 ^ P(hj ^)j(xo ' z o) II P 2 ^ h J x k^yko - z k 0 )^ 


/2 


or 


i-1 


" P(hj xp)I Czo ~ £. 0 ) 
0=o 


i-1 


■1“ * o I P (k J x k ) 1)1 y ko " 2 ko I (37) 


depending on which norm is used. 


* • 

Now ReX k < o, k=1,...,n; thus, there exists h k > o such that | P(^k^k ) I = 1 » 

« 

for each k and for all h in the closed interval o < h < h^. It follows that 
|p(X k h)| < 1 


Set 


* 

h* = min , (h* is positive) 

k 


(38) 


N-1 

Then if the sequence {h j} ^ has the property that hje(o,h*), for all j, 

it is concluded that |p(hjX k )| < 1 for all j and k, and consequently, 


n /i-1 \ \ 1/2 

l v p 2 (hjX k ) (yko - z ko> 2 ) < II Zo - z o 

^k=1 \j=o 


and 


12 
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T \/ 0 ' P(hjXk> l y ko “ 2 ko I < II lo - lo II 


Equations (37) then give the necessary conclusion of stability of the numerical 
equation (24). 

Since amplification factors and matrices have dominant roles in the study of nu- 
merical stability, it is of value to examine their graphs in detail. 


4.1 AMPLIFICATION FACTOR GEOMETRY 

Graphs of the region of stability (e.g., fig. ^) are too coarse for the present 
analysis. This section will consider only a four-stage, fourth-order explicit 
Runge-Kutta method. 

In this case, 


4 

p(hX) = Z 

J=° 


(hl)J 


where X is complex. 


In polar coordinates, 


X = p(cos 0 + i sin 0) 


and 


4 (hp(cos 0 + i sin 0)1J 
p(hX) = Z 

jso J* 

or 


p(hX) 


4 (hp) J . , 

Z (cos(j0) + i sin(j0)J 

j=o 


(39) 


(40) 


where equation (40) resulted from applying De Moivre's theorem to equation 
(39). Thus, 


13 
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( p(hX ) | = 



(hp)J 

“JT 


C03 ( J0 ) 



(hp)J 

~7T 


ain(je) 


2 


1/2 


(41) 


The first point of view will be to specify various choices of X and plot 
|p(hX)| as a function of h, h > o. Figures 3 through 7 each consists of three 
plots of |p(hX)| with the angle 6 fixed, but with p = 2,5, and 10. Thus, 
each figure represents the evolution of the graph of |p(hX)| as the modulus 
p of X Increases from 2 to 10 (for a fixed angle). However, by fixing p 
(e.g., p=2) and selecting the corresponding plot in each figure, it can be seen 
that evolution of [p(hX)| as 0 decreases through the range of 180° (fig. 3), 
120° (fig. 4), 90° (fig. 5), and 85° (fig. 7). Figure 6 (X=o) is included 
because it clearly represents a limiting case of the purely imaginary situation 
(fig. 5); i.e., 


| p( o ) | = lim |p(iph)| 
p o 


This indicates that, from the numerical stability view, the behavior of problems 
with X=o (as in ref. 4) is more closely related to problems with small, purely 
imaginary eigenvalues than with small, real negative eigenvalues. The concave- 

downward segment of the graphs for the purely imaginary case (0 = 90°) suggests 

that small perturbations could produce either stability or instability if the 

step sequence were confined to certain intervals (i.e., if X=2i, then 

hi e (0,0.5)). 

Additional global Information about the stability region S can be obtained by 
setting p = hX in equation (41) and finding contour lines 


] p(p) i = r, o < r < 1 


Figure 8 illustrates the level curve configuration for the fourth-order Runge 
Kutta methods, with r=0.1 k, k=1,2,...,9. To obtain both quadrants, reflect 
figure 8 across the negative real axis. 


4.2 NUMERICAL KERNELS 

While theorem 1 applies only to diagonal systems (eq. (29)), it is possible to 
bring more general systems under that theorem’s sphere of influence. Consider 
a system 


i. = Ax 


(42) 


14 
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where A is an n by n matrix and such that A is diagonalizable . Thus, there 
exists a nonsingular matrix M, and 


M" 1 AM = 


( 43 ) 


n 

where ^ is diagonal. If the eigenvalues of A are {X^}^ ^ , then it is 
well known that 


% = diag {X , X n } 


( 44 ) 


(i.e., eigenvalues are invariant under a similarity transformation). 

Remark. The order of the eigenvalues in is dependent on 4 . 

Amplification matrices and numerical kernels for equation ( 42 ) is now considered. 
In fact, the development in equations ( 1 ) through (8) extends naturally to vector 
systems. Thus, if 


f(x) = Ax 


( 45 ) 


then equations (6) become 


where 


£1 = 

A*n 



£2 * 

a(i 

+ a 2 (h n A) ) x n 


£3 = 

a(i 

♦ a 3 (h n A) + a 2 b32 (h n A) 2 ) x^ 


£4 = 

a[i 

+ a4 (h n A) + (a 2 bi|2 + 23843) (h n A) + 2^32643 

(h n A) ^ ) Xjj 

is 

the 

n by n identity matrix and 



= P(h n A)x n 

( 46 ) 


where 


15 
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P(h n A) s I + h n A ♦ ^ (h A)* ♦ (a 6 0)3 + 82^42^4 

* 83643414) (h n A)^ + (82632643414) (h n A)** 


(47) 


Farther, to 0(tr), 
n 


h n A 


(h n A)^ 


j=o 


J! 


which implies that the matrix analog of equation 03) is namely 


p(h n *) = £ Lisiii 

j=o J* 


(-48) 


N-1 


If the view is taken that a step-size sequence { h 4 } is given, then at any 

i=o 

i-1 


time tj_, tj_ * t 0 + £ hj, the i th iterate is defined recursively by 

j=o 

formula equation (46) or 


x L = P(h i _ 1 A) P(h i-2 A)...P(h 0 A) x< 


(49) 


or 


i-1 

IT 

j=o 


x, =( IT P(hjA) 1 Xq 


( 50 ) 


where it is understood that the it product in equation (49) is ordered as in 
equation (49). The amplification matrices in equation (49) are not commutative. 
The it product in equation (50) is called the numerical kernel (associated with 
the i th iterate). An important relationship is now established. 

Observation 1 . Amplification factors (for Runge-Kutta methods) are invariant 
under similarity transformation (i.e., if A is similar to then P(h n A) 
is similar to P(h n 4>n)). 


16 
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Proof. First , 


A 1 = AM(M _1 AM) 1 " 2 M" 1 A, for i > 2 


(51) 


but 


4 (h A) 1 

P(h n A) = z -2-— (52) 

i=o i' 

and consequently, 

4 hj, 1 (M" 1 AM) 1 

M" 1 P(h n A)M = Z ~ 

i=o 11 

t ( h n W 1 

2 1=0 li 

= P(h n 1|> M ) (53) 


This completes the proof. 


It should be observed that the subscript of is necessary in the sense that 

while the diagonal elements of are exactly {A A n ), the order may not 

be preserved; i.e., 4 *m = diag t^il »^i2» * • ' »^in^ where the two eigenvalue sets 
are equal in elements and multiplicities. 

The previous result is extented to numerical kernels. 

Observation 2 . Numerical kernels for Runge-Kutta methods are invariant under 
similarity transformation (i.e., if A is similar to t then K (A) is 
similar to K i (^ M )). 

Proof. From equation (50), the numerical kernel for equation (42) is 


i-1 

K a (A) = it P(hjA) (54) 

j=o 
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Hence , 


M" 1 Ki(A)M = 


i-1 

n 

J=o 


(hHPthjAjM) 


i-1 

= n P(hj4»M) « Ki(^ M ) (55) 

j=o 


Observe that the numerical kernel for is diagonal. Equation (54) introduces 
the notation for the kernel; thus, the i th iterate of a numerical oolu- 

tion starting at with step sequence Ih^; is 

i=o 


Xj - K^( Ajx^ , i=1 ,2,.. . ,N 


(56) 


Equation (56) generates a numerical approximation to 


x = Ax , x (t 0 ) = Xq 

By multiplying equation (56) by M -1 on the right, it follows that 


M-’xi = M" 1 Ki(A) x 0 

= M" 1 Ki(A) MOT 1 x D ) 


M _1 3Ci = Ck ± OP m )) (M- 1 x 0 ) (57) 


where observation 2 was used to obtain equation (57). 
If 


lo = M_1 *o 


is set in equation (57), then that equation generates an approximate solution 
to the initial value problem 
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i = ^mx. » y<t 0 ) * &> 

Multiplication of equation (57) by M (on the left) gives 

ii = (M K^m) M" 1 ) x 0 (58) 

It should be noted that K^('I'm) and Kj_(A) are step sequence dependent; i.e., 

. sN-1 

functions of ihu 
4.3 STABILITY BOUNDS 

Theorem 2 . Consider the initial value problem 

£ = , y. (t G ) = i 0 ’ (59) 

where ^ « diag {X lf ...,X n } with Re Xj < o for all i. Let 



be a numerical solution of equation (59) generated by an explicit Runge-Kutta 
algorithm. Then there exists symbols P^, Tj, and an interval (hL,hRj such 
that o < P k 1 1 and o < < 1 , and if hje(hL,»hRj for all i, then the 

k th component of £i, namely y^, satisfies the inequality. 


T k I Yko “ z ko I - I Yki “ z ki I 1 ^k I Yko “ z ko 


for all i , 



1 = 1 , 2 ,..., 



Proof. For each 


N and all k, k=1,...,n. The inequalities in equation (60), 
is any other numerical solution. 

there exists h^ > o such that 


I P( x k h k ) | = 1 
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and for all h, o ^ h < h k 

I P(*k h > I < 1 


Set 


h # = ain {fy} 
k 


and suppose 

(h L ,h R ) C (o,V) 


where 


hL = min {h*}, h R = max {hi) 
i i 


Since 


fjj(h) = | p(Xifh) | 


is a continuous function and (hL,h R ) is a compact set, there exists 
symbols xic and such that 


min | P(Xkh) | = Tk < | p(Xk^) | < p* = max | p(Xkh) I 
he(h L ,h R J he(h L ,h R J 

where the inequalities in equation (61) hold for all he(hL,h R ). 
Now 


/i - 1 

li = " 
\J=o 


p < h j40) Zo 


positive 

( 61 ) 
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which implies that 


fi-1 


y ki P^k h jy y ko 


(62) 


and 


i-1 

z ki a * PU k hj) z ko 
J*o 


(63) 


Subtracting equation (63) from equation (62), and taking absolute values 


fi-1 

y ki - z ki l a l * | P<Xkhj) (][ y ko * z ko ! 

J*o 


(64) 


but 


or 


i-1 i-1 i-1 

n T k < it | p(X k hj) | < it p k 

j=o j=o J*o 


i i-1 i 

T k < TT | p(X k hj) | < p k 
j=o 


which, with equation (64) implies that the conclusion for each i and k. 

Discussion . The interval (hL,hp) is referred to as the step interval for a 
particular numerical solution. Since J p(X k h) | * J p(X k h) |, then if conju- 
gate pairs of complex eigenvalues occur, it is only necessary to compute r k 
and p k for one of the pair. Figure 9 illustrates the theorem for equation 
(59) with eigenvalues { +5i, -2, 5(cos 120° + i sin 120°)). Clearly, the bounds 
would be sharper if step variations were small; that is, hR - h^ were near 
zero. 
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5.0 ERROR PROPAGATION FOR HONDIAGOHAL SYSTEMS 


From the study of Lyapunov stability (ref. 5), It is well known that if a 
linear, constant coefficient differential system 

x * Ax 

has eigenvalues {(Xi)i«i) and Re Xi <0 for all 1, then the solutions are 
stable (in fact, asymptotically stable), «nd this is the case even if Mere are 
multiple eigenvalues. An analogous question in nunerioal stability would be the 
following: Consider two constant coefficient, linear systems of equal dimension 


x = Aix 


£ ! *2J 


in which the Jordan canonical form of A^ and A£ are given by: 
Mi 1 A! Mt s J v , Viz *2 * J 2 


and the eigenvalues of and J 2 are identical in form and multiplicity 
but is diagonal and J£ is not. Thus, 


* diag (x^ , • . . ,X^ ! X2 > * • * 1X3 i • • • »X|» , • • • ,Xp) 


where 


r 

Z multiplicities of Xi = n, 
i = 1 


and 


J 2 = diag (A 1t A 2 , — ,A r } 


where A< is a square matrix with dimension equal to the multiplicity of 
and (ref. 6) 
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6 iJ 


o or 1 


Then, as has been proven in theorem 1 , if Re Xi < 0, the numerical solutions 
of 


X = 2 lX 


-1 N-1 

where jr = M] x are stable provided the step sequence ^ 


is admissible in 


the sense of the statement of the theorem. What then can be said for the stabil- 
ity of the solutions of 


X = J 2 X 


Since the test algorithm is fourth order Runge-Kutta, a related question con- 
cerns the relationship, if any, between the contour curves for |p(y)| (fig. 8 ) 
and error propagation (i.e., does the topography of |p(p)| influence error 
propagation?), figure 8 has several distinct features including two wells 
and a "plateau" centered on the imaginary axis and extending out about 1.5 units. 

Specifically considered are examples having the general form 


J 2 


Xi 

6 1 

0 

0 

o 

Xi 

62 

0 

0 

0 

Xi 


o 

0 

0 

Xi 


(65) 


where Xi is the complex conjugate of Xi- Equation (65) can include complex 
eigenvalues of multiplicity two (e.g., Xi = i, Xi = -i where i 2 =- 1 ) and 
nonpositive real roots of multiplicity four. These cases are of particular interest 
in orbital motion equations and mechanical systems. The amplification matrix 
associated with equation ( 65 ) is 


o 



4 

P (h J 2 ) = 2 

J*o Jl 


and the error vector e» propagates according to 



Sm+1 = p <W 2 ) ha 


79FM32 

( 66 ) 


where 


ha = ( e 1» 


e m ) ^ 
• • » c n ; 


(67) 


and a constant step size h is assumed. The amplification matrix P(hJ 2 ) 
is given by 


P (h J 2 ) =1 


^p(hXi) Pi2 Pi 3 Pin 
o p(hXi) P 23 P24 

o o p(hX’-j) P34 

\.o 00 p(hX^ 


P12 


P23 


( 68 ) 


where p(hX) is defined in equation ( 13 ), 

= h fl + h Xi + ( h X 1 )2 + ( h x t )3l 
V 2! 3! ’ 

- h 6 2 P ■*" h(X-| + ) + h (Xi + X 1 X 1 + Xi ) + 

V 2! 3! 

h 3 (X| + x 2 + iiX 2 + X^ ) \ 

4! / 

P34 = h 63 + h X-j + (hXi) 2 + ( hX i ) ^ ) 

V 2! 3 ! / 

p 1 3 = Y? 6 -|< 5 2 + h (2X1 + X-| ) + h^( 3 X^ + 2 X^ + X^) j 

\2! 3i 4» / 


and 
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P14 


= 6 


1* S 2 <5 3 + h . (fftl -t 


In the study of error propagation, the equations to be considered are: 


m +1 


e“"' = p( hX t ) e“ + p 12 e 2 + P 1 3 £3 + Pi4 £4 


m 


m 


e 2 +1 = p(hX^) e 2 + P 23 £3 ♦ p 2 4 £? 


m+1 

£ 3 


= p( hX -j ) e“ + P34 eij 


e^ + 1 = pOiX^) ejf 


(69) 


It will be assumed that control of error propagation has occurred if for any 
vector e° = (£°, £°, £°, e 0 )^ with components le°| < 1 , it follows that the 
1 2 3 H i ~ 

components e ” +1 of e 0 * 1 (defined iteratively by equation (69)) satisfy 


tef 1 ! i 1 


(70) 


for all i and m. 

This can be phrased as: all components of the error vector iterates are 

within the unit disk if the initial error has components within that disk. 

In this document, the analysis is restricted to equation (65) with 61 =1 , 6 2 
and 63 = o. An examination of pjj factors shows that p^j = o except 
for i = 1 and j = 2 and in this case, 


PI 


/ ( hX 1 )2 (hX,)3\ 

2 = h 1 1 + hX 1 + + ) 

V 2 ! 3'- / 


( liX t + (hXj)^ + ( hX 1 ) 3 + ( hX 1 ) , if Xi | o (71) 

21 3 ! 


Note that 
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lim |p 12 | = h 
|Xi| + o 

That is, the graph of (pi?! (as a function of h) approaches the 45° line 
as JXtI >o. 

Equations (69) reduce to 

m+1 ... v m m 

ei = p(hX-| ) ei + P 12 £2 

£j +1 = p(hX) ej , j = 2, 3, 4 and X = Xi or (72) 

In contrast to our previous study of amplification factors for diagonal matrices, 

the factor e° +1 is propagated as coupled function of e? and ef, the amplifi- 
cation factor P 12 is a function of the eigenvalue X and finally, the triangle 
inequality is too coarse to use in the analysis; i.e., 


|e? +1 | = | P(hX-|) e? + P12 e? | 

< | p(hX;) | | e“ | + | P12 | | ef | 


does not admit the type of conclusions desired. Certain conclusions do seem jus- 
tifiable for eigenvalues of the form i 3 , 3 :> o from an analysis of | p 12 1 » 
especially when coupled with the previous work on |p(hX)|. Figures 10 and 11 
give the plots of |pi 2 | for X= p(cos 0 + I sin 0 ) , with p= 2 , and for selected 
values of 0. When the graph of jpi2| for 9 = 90° (i.e., Xi=2i) in figure 10 
is compared with the graph of |p(Xh)| in figure 3, it i 3 apparent that error 
propagation will not be controlled. This observation also holds for X=o; 
e.g., zero is eigenvalue of multiplicity four but with nondiagonal Jordan form. 


It should be realized that it is unnecessary to be concerned with gj , j > 2 
since we have assumed that hX is within the stability region S, and conse- 
quently 
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Preliminary work has indicated that for certain initial errors (£j = 1, j 
s 1,2 ,3 ,**) there are regions in the second quadrant for which the first iteration 
fails to satisfy inequalities in equation (70); i.e., 



Figure 12 illustrates the situation for X = 2(cos© + i sin©), and the polar 
coordinate distance is Ph = 2h. 


A comparison of figures 8 and 12 suggests that the plateau centered on the imagi- 
nary axis had significant influence on the plot of | j . It has also been 


observed that if iterations of £° are continued according to equation (72) 
then equation (70) may be met for sufficiently large m, even if operating in 
the area between the imaginary axis and the boundary of | e 1 1 < 1 . 


6.0 CONCLUSIONS 


The plots of stability regions alone are inadequate to study stability or error 
propagation of numerical solutions that correspond to linear systems having re- 
peated complex roots, and the problem is more acute for eigenvalues near the 
imaginary axis; i.e., with angles slightly larger than 90°. There is increasing 
evidence that most numerical integrators such as explicit Runge-Kutta and linear 
multistep methods have difficulty with such differential systems and improvement 
of the performance of the algorithms will depend on the recognition of what 
causes the problem. It is suggested in this paper that more emphasis must be 
placed on error propagation in the case that Jordan Canonical form of the A 
matrix is nondiagonal. That is, where the system is x = Ax. In fact, a lack 
of understanding of exactly what numerical stability is or should be is a large 
part of the problem. In this paper, various definitions are put forth for numer- 
ical stability, and the amplification matrix is studied in detail (for the 
fourth-order Runge-Kutta method) in order to understand reasons for the apparent 
failure of stability in the situations given in this paper. For the case in 
which the A matrix is diagonalizable, a fairly complete theory exists and sev- 
eral theorems are proven for that case. Theorems on the existence of numerical 
stability and stability bounds are also presented. 

Evidence that a plateau, centered on the imaginary axis and extending out from 
the origin, may have sufficient influence to cause instability for repeated 
eigenvalues when arguments 6, 90° < 0 < 120° are given. It is also noted that 
(from the numerical stability viewpoint) an eigenvalue of zero should be con- 
sidered as a limiting case of i8 where i^ = -1 and B + o rather than as a 
negative real number (tending to zero). It is clear that additional studies are 
required with computer support, but the preliminary analyses by the techniques 
of this paper appear to have potential. 
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Figure 1.- Region of absolute stability of fourth-ordei explicit 
Runge-Kutta method (ref. 7,p.41). 



^0* A^y/ 

radius *"0 = 5 


Figure 2.- A "discrete tube" t x R 2 associated with a numerical solution. 
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Stability zone 
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Figure 5.- Modulus of amplification factor versus step size (x = pi, for p = 

2, 5, and 10 and i 2 = -1). 


Stability zone 
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Figure 7.- Modulus of amplification factor versus step site (i = ,,{cos 85* 4 

i sin 85°) for p = 2, 5, and 10). 


Ui 

cn 



Figure 8.- Contour curves for | p( y ) J where p{ y ) is 
for fourth-order Runge-Kutta method. 


the amplification factor 


|p(*h)J 
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Figure 11.- Modulus of amplification factor p, 2 versus step size for 
> - 2(cos e + 1 sin e) where e ■ 140”, 160°, and 180°. 
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